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Abstract 

The parametric instability arising when ordinary differential equations 
(ODEs) are numerically integrated with Runge-Kutta-Nystrom (RKN) 
methods with varying step sizes is investigated. Perturbation methods 
are used to quantify the critical step sizes associated with parametric 
instability. It is shown that there is no parametric instability for linear 
constant coefficient ODEs integrated with RKN methods that are based 
on A-stable Runge-Kutta methods, because the solution is nonincreasing 
in some norm for all positive step sizes, constant or varying. 

1 Introduction 

In a recent paper |10j . Wright showed how instability can arise when linear 
second-order constant-coefficient ordinary differential equations (ODEs) are nu- 
merically integrated with varying step sizes, even when all steps are smaller than 
the critical step size that ensures stability in constant-step computations. One 
of his examples was the central difference method applied to the model problem 

x + x = (1) 

With constant step size h, the central difference method is known to be stable 
when h < 2. Wright showed how it can become unstable when the step size 
has small-amplitude oscillation of period 2 about the constant value h = \/2. 
The instability manifests itself as oscillation with growing amplitude. Instability 
also arises if the step size has period-3 oscillation about h = 1, and so on to 
period-6 oscillation about h sa 0.518. He conjectured that instability could arise 
for arbitrarily small step size. Shortly after Wright's paper was published, Skeel 
wrote a letter to the editor [5] pointing out that he had already proved Wright's 
conjecture in a short note published in BIT [TJ. 
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In this work the instability phenomenon is investigated by two approaches: 
first, stability analysis using contractivity concepts, and second, parametric in- 
stability analysis using perturbation methods. Here a brief description of these 
approaches is given. 

In the usual stability analysis of a numerical ODE integration method, the 
linear ODE system is transformed to diagonal form using a similarity transfor- 
mation. Because the transformation depends on the the step size, this approach 
is only useful for methods used with constant step sizes. Contractivity analysis 
is not based on the diagonalising transformation. The idea is, given an ODE 
problem whose solution is nonincreasing in some norm, to derive conditions un- 
der which the numerical solution will also be nonincreasing. Hairer et al. [2 
have used contractivity concepts to derive stability results for a wide class of 
ODE integration methods with constant step size. The same approach is used 
here to derive stability results for methods with varying step sizes. It is shown 
that A-stable methods are stable also with variable step size. 

The second line of inquiry is prompted by Wright's observation that 

. . . integrating a variable-stiffness system with constant time steps 
resembles integrating a constant-stiffness system with varying time 
steps. 

It is well known that oscillatory variation of the stiffness parameter in variable- 
stiffness systems can induce a kind of instability called parametric resonance [5] . 
Wright's remark suggests that the instability he observed is a kind of parametric 
resonance of the "discrete time system" comprised by the numerical integration 
method applied to the ODE. Indeed, Wright's stability region diagrams resemble 
the classic Strutt diagrams used in analysis of parametric resonance. 

Although the literature on parametric resonance is vast, there appears to be 
little on parametric resonance of difference equations. Tanaka and Sato [5] have 
studied a second order difference equation with periodic coefficients, a kind of 
discrete Mathieu equation. They compute Strutt diagrams and show how to 
approximate the stability region boundaries using perturbation methods. This 
type of analysis is applied here to RKN integrators with periodically varying 
step size. 

2 Runge-Kutta-Nystrom methods 

In this section a well-known general class of numerical integration methods for 
second-order ODE initial value problems is presented. Five specific methods 
from this class are singled out for further study. 

The s-stage Runge-Kutta-Nystrom (RKN) method for advancing the solu- 
tion of the second-order velocity-independent ODE 

x = f(t,x) (2) 

from t = t n to t n+ \ = t n + h n is 

s 

ki = f(t n + Cih n , x n + Cih n x n + h n S ' ciijkj) (1 < i < s) 
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For example, the coefficients for an explicit fourth-order RKN method [J, p. 262] 
are 
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A RKN method is equivalent to a Runge-Kutta method if there exists a 



coefficient matrix a,j such that 
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For example, the A-stable third-order 2-stage SDIRK method [31 p. 203] corre- 
sponds to the RKN method 



(4) 



where a — (3 + \/3)/6. 

The Newmark method for @ is a RKN method with two coefficients, (3 and 
7, that appear in the RKN tableau as follows: 
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The central difference method is given by the Newmark parameters 

P = 0, 7 = \ 



(5) 



The Newmark method that was shown by Wright to exhibit parametric insta- 
bility is specified by the parameter values 



P = l 



1 



(6) 



The only choice of Newmark parameters that gives a Runge-Kutta method (the 
trapezoid rule) is 

P=\>1=\ W 



1 In [6] it is shown that no R-stable third-order diagonally implicit RKN method can have 
0.547 < ci < 1.213, but this is only for RKN methods satisfying the simplifying condition 
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Applying the RKN method to the linear constant coefficient model prob- 
lem ((IJ and eliminating the stage variables gives 



y n+1 = R(h n )y n 
where the state vector y n is defined as 



(8) 



tin 



and the one-step transition (amplification) matrix R is given by 
R{h n ) -- 



1 - h*F(I + hlAY^e h n - hlb T (I + h 2 n A)- x c 
-h n b T (I + hlA)~ l e 1 - h 2 n b T (I + hlA)- l c 



with e being a vector of ones. In particular, the transition matrices for the RKN 
methods presented earlier are as follows. 

• The explicit fourth-order RKN method given in Eqn ([3]): 

-l/2 + ft 2 /24 -ft/6 



1 ft ' 




-ft 1 


+ h 2 



R(h) , , 

- ft/6-ft 3 /96 -l/2 + ft 2 /24 

The A-stable third-order 2-stage SDIRK of Eqn (g}: 
R(h) = 



1 ft 
-ft 1 



+ - 



(1 + h 2 a 2 ) 2 
The Newmark method 
R(h) = 



-(4ft 2 a 3 -ft 2 a 2 + l)/2 ~(h 2 a 3 - a + l)ah 
{h 2 a 3 -a + l)ah (Ah 2 a 3 - h 2 a 2 + l)/2 



1 ft 
-ft 1 



-1/2 -/3ft 
7/1/2 -7 



(l + ft 2 /3) 

with, as special cases, the central difference method of Eqn (JSJ: 

-ft/2 
ft/4 -1/2 
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R(h) = 

the Newmark method of Eqn © 
i?(ft) 



1 ft 
-ft 1 



(l + ft 2 /2) 
and the trapezoid method of Eqn (0 : 

ft 2 



i?(ft) = 



1 ft 
-/i 1 
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3 Contractivity Analysis 



The main result of this section is the following theorem, which states that A- 
stable Runge-Kutta methods such as the trapezoid rule are contractive in a 
norm that is independent of step size. The same result of course holds for RKN 
methods derived from A-stable Runge-Kutta methods. Thus, these methods do 
not exhibit the parametric instability caused by varying step sizes. 

Theorem 1 For any A-stable Runge-Kutta method applied to a stable linear 
constant- coefficient ODE there exists a symmetric positive definite matrix W , 
independent of the step size, such that the numerical solution is contractive in 
the W -weighted euclidean norm, that is, any two numerical solutions ymVn 
satisfy 

WVn+l - Vn+lWw < \\Vn ~ Vn\\w 

Before looking at the proof of the theorem, a couple of remarks. First, 
Wright's results for the implicit R-stable Newmark method of Eqn ^ is a 
counterexample showing that R-stability or implicitness are not sufficient to 
preclude parametric resonance with varying step size. 

Secondly, in his letter to the editor, Skeel [8] mentions that the Newmark 
method of Eqn © 

. . . does not satisfactorily deal with variable step size. A method that 
does is the implicit midpoint method, for which the amplification 
matrix is an orthogonal matrix. Of course, it has the same practical 
drawback of implicitness. . . 

The implicit midpoint method has the same transition matrix as the trapezoid 
rule of Eqn |7p. Because the transition matrix is orthogonal, the method is 
contractive in the euclidean norm. Thus, theorem 1 extends Skeel's remark to 
cover all RKN methods that are equivalent to A-stable Runge-Kutta methods. 
However, since A-stable methods are implicit, this result is of limited use to 
those who want to use explicit methods. 

The theorem is proved using three lemmas. The first lemma states that any 
stable linear constant-coefficient system is contractive in some norm. 

Lemma 1 If A is Hurwitz then there exists a constant symmetric positive def- 
inite matrix W such that the solution of the ODE y = Ay is contractive in the 
W -weighted euclidean norm. 

Proof. R suffices to consider the homogeneous problem 

V = Ay 

Since A is Hurwitz, there exists a symmetric positive definite matrix W that 
satisfies the Lyapunov equation 

A T W + WA = -I 

Apply the linear transformation z — W x / 2 y to take the original ODE into 

i = W^AW-^z 
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The solutions to the ODE are monotonically nonincreasing in the euclidean 
norm, because 

^-Ml 2 = -.^w- 1 * < o 

at 

Since \\z\\ = \\y\\w> the stated result follows. ■ 

The s-stage Runge-Kutta method for advancing the solution of the linear 
constant coefficient ODE 

V = Ay 

from t n to t n+ \ = t n + ft, is given by 

s 

9i = Vn + h^ciijAgj (1 < i < s) 

S 

Vn+i = y„ + h^bjAgj 

Eliminating the stage variables gi results in 

y n+ i = R(hA)y n 

where R is the stability function of the Runge-Kutta method. 

A linear transformation z = Vy, where V is a constant nonsingular matrix, 
gives a new ODE 

z = VAV~ x z 

The next result tells us that the result of applying the Runge-Kutta method to 
the transformed ODE is the same as transforming the numerical results from 
the original ODE. 

Lemma 2 The following diagram commutes: 





V 




Vn 




Zn 



R{hA) i i RihVAV' 1 ) 





V ^ 




y n +i 




Zn+l 



PROOF. Substitution. ■ 

Finally, the following result from [1J Corollary 11.3] is needed. 

Lemma 3 // an A-stable Runge-Kutta method is applied to a linear constant- 
coefficient ODE that is contractive in the euclidean norm, then the numerical 
solution is contractive in the euclidean norm. 

Combining the three lemmas gives the theorem. Lemma 1 gives the weight 
W under which the ODE is contractive. V = W 1 ^ 2 is the linear transformation 
matrix that relates the original variables y to the new variables z that are con- 
tractive in the euclidean norm. By lemma 3, numerical sequences z n computed 
with A-stable Runge-Kutta methods inherit the contractivity of the ODE, and 
by lemma 2, this also applies to the original ODE. 
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4 Stability charts 



A necessary and sufficient condition for stabie integration of the model prob- 
lem ([!} with constant step size is given by the Schur-Cohn inequalities 

|trace(i?)| - 1 < det(R) < 1 (9) 

Thus, with constant step sizes the central difference method, Eqn ([5]), is stable 
for < h < 2, the RKN method of Eqn © is stable for < h < 2(2 + 2 1 / 3 - 
2 2/3)i/2 _ 2.59, and the SDIRK method of Eqn flU, the trapezoid method 
(Eqn ((?])), and the Newmark method of Eqn © are stable for all h > 0, i.e. 
they are R-stable. 

Now consider integration of the model problem ([1]) with a sequence of step 
sizes h n that is periodic with period p. The composition of p steps can be 
interpreted as a one-step method integrating from t = t n to t = t n+p = t n + 
h n H h /i„+ p _i. The application of one step of the composed method gives 

Un+p = R{[h n , ■ ■ ■ , h n + p -i])y n 

where the transition matrix R of the composed method is given by the matrix 
product 

R([h n , . . . , h n+ p-i]) := R(h n+p _i)R{h n+p -2) ■ ■ ■ R(h n+ i)R(h n ) 

This composed transition matrix determines the stability of the linear homoge- 
neous difference equation ©. The stability condition is given by the inequali- 
ties ^) applied with the composed transition matrix. 
For the periodic step size variation given by 

h n = h + ecos(27m/p) (10) 

a stability chart can be made by numerically testing the stability criterion for 
a large number of specific values of (h,e). Figure 1 shows values giving an 
unstable composed Nystrom method of Equation when the step size period 
is p = 6. Also shown are the corresponding stability charts for the central 
difference method and the Newmark method of Equation©. 

The instability regions for the central difference and Newmark methods form 
wedges with their points on the /i-axis. Stability charts for other periods p are 
similar, except that the number of wedges increases with p. The figure resembles 
a Strutt diagram for the Mathieu equation [5]. 

The instability regions for the Nystrom method are more rounded and, for 
values of h smaller than the constant-step stability limit, do not reach the h-axis. 
This figure resembles Strutt diagrams for damped oscillators. 

The stability charts for the trapezoid and SDIRK methods are not shown, 
because they contain no instability points. This is because they are A-stable 
Runge-Kutta methods, which according to the theorem in section [3] remain 
stable with varying step size. 

5 Perturbation Analysis of Stability Regions 

The points on the stability charts where the instability regions intersect the h 
axis will be called critical step sizes. In this section, their values are identified 
using perturbation analysis. 
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Substituting the trial solution 



Vn 



and the step size (|10l) into © and equating powers of e gives a sequence of 
coupled difference equations, the first three of which are 

y^i = R(h)yW (12) 
ifil = R(h)yP +co S (2nn/p)R'(h)yW (13) 
/£> = R{h)yW + cos(2t™/p) + W(h)y^ 



Vn+l -~v/i*n 1 — ■-/ x-y i -~ v/an 1 ^' 

The first equation, Equation (|12p. is a homogeneous linear constant coefficient 
difference equation, and is stable provided that the constant component h of 
the step size is within the stability range for the constant-step method. 

The remaining equations are linear constant coefficient difference equations 
with a harmonic forcing term. If the eigenvalues of R are smaller than one in 
magnitude (i.e. the method has algorithmic damping), then these equations will 
all have bounded solutions. For such methods, the "wedges" in the stability 
chart do not reach the h axis, and the method is stable when the amplitude e of 
step size oscillation is sufficiently small. This is confirmed by the stability chart 
(Figure 1) for the Nystrom method of Eqn ([3]), which has algorithmic damping. 

Resonance instability may arise in integration methods without algorithmic 
damping if the forcing frequency matches the natural frequency. It therefore 
suffices to restrict attention to methods whose transition matrix R(h) has eigen- 
values of unit modulus. These eigenvalues can be written as 

eig{R{h)) = exp(±iw(h)) (14) 

where w is a real function of h. For example, the central difference method has 



cj(h) = arctan(/i-v/± - /i 2 /4, 1 - h 2 /2) (0 < h < 2) 
while the Newmark method of Eqn © has 



uj(h) = arctan(/iv / l + h 2 ) (0 < h) 

For methods satisfying (|14l) . the solution of (fl"2"|) is a linear combination of 
trigonometric functions of the form 

= Yi cos(nw) + Y 2 sin(no;) (15) 

Substituting (|T5]) into (fT3"|) gives a forcing function for y„ of the form 

R'(h)Yi cos(nw) cos(27rn/p) + R'(h)Y2 sin(nw) cos(27m/p) 

Resonance occurs when 

u(h) = - (16) 
P 

or 

u){h) = tt 

P 
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p 


h° 


/i 1 


2 




1.4142 




±1/2 




3 


1 






±1/8 




4 




d 0.7654 


±(2 


-V2)/8R 


i ±0.0732 


5 


(V5-1)/2p 


s 0.6180 


±(3- 


V5)/16r 


i ±0.0477 


6 


(v/6- v^)/2p 


d 0.5176 


±(4- 


V / 12)/16k 


i ±0.0335 



Table 1: Smallest critical step size h° and wedge width h 1 for period-p oscilla- 
tion, central difference method 
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h° 
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00 








3 




d 1.4142 


±1/4 




4 


V-2 + 2V2R 


d 0.9102 


±(-1 + V2)/4p 


a 0.1036 


5 


V-4 + 2V5R 


d 0.6871 


±(-2 + V5)/4s 


d 0.0590 


6 


|V-18+12V3r 


i 0.5562 


±(-3 + 2V3)/12? 


d 0.0387 



Table 2: Smallest critical step size ft and wedge width /i 1 for period-p oscilla- 
tion, Newmark method 



Values of h satisfying these equations are critical step sizes. 

The smallest critical step size for period-p step size oscillation is denoted h°, 
and is given by the solution of Equation (1161) . The minimum critical step sizes 
for various p are listed in Tables [T] and [2 These values agree with those found 
by Wright [ID] by other means. 

Any RKN method of order at least 1 has 



R(h) 

and so 



1 h 

-h 1 



0(h 2 ) 



eig(R{h)) =±ih + 0(h 2 ) 
The smallest critical step size for large p is then given by 

min h CTit a ^ 

Thus, the critical step size can be arbitrarily small. This agrees with Skeel's 
result [7j for the central difference method. 

The method of strained parameters is now used to analyse the shape of the 
instability region "wedge" near the /i-axis. The average step size parameter is 
assumed to have the form 

h = h° + h 1 e + 0{e 2 ) 
where h is a critical step size. The step size then varies according to 

h n = h° + h x e + e cos(27m/p) + 0(e 2 ) 
Substituting this into the composed transition matrix gives 

]) - R + eRt + 0(e 2 ) 
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where 

Ro := R(h ) 3 

and 

Ri := (cos{2nO/p)+h 1 )R(h )P- 1 R'(h°)+. . .+(cos(2 7 r(p-l)/p)+/i 1 )i?'(/i )i?(/i )^ 1 

The boundaries of the instability regions are found by solving the stability con- 
ditions for the composed transition matrix. 

For example, for the central difference method with step size period p = 3 
we have h = 1 and 

R([h 2 MM) = - / + ^ 

Neglecting terms in e 2 , the stability condition ^ is 

l<l-^e 2 + 12e 2 (^) 2 <l 

which is satisfied when 

8 

The instability region boundary is therefore approximated by 

h=l± l -e + 0{e 2 ) (17) 

Figure 2 shows instability points in the neighbourhood of the critical step size 
for the central difference method for parametric oscillation of period p = 3, and 
the instability region boundaries predicted by (I17|) . The approximation appears 
to be correct. Tables [1] and give instability region "wedge" widths for the 
central difference method and the Newmark method. 



6 -Mh 1 - 4 
48/j 1 - 3 -6 



+ 0(e 2 ) 



6 Conclusions 

The following results were presented: 

• A-stable Runge-Kutta methods are stable when applied with varying time 
step to stable linear time invariant ODE systems. 

• For numerical ODE integration methods that have numerical damping, 
there will be no parametric resonance if the step size oscillation amplitude 
is sufficiently small. 

• The critical average step size about which small oscillation may cause 
parametric resonance in the numerical integration of a simple oscillator 
can be found using perturbation methods. 

• The perturbation analysis indicates that there can be arbitrarily small 
critical step sizes. 

• The method of strained parameters can be used to approximate the shape 
of the instability region in Strutt diagrams for small amplitude step size 
oscillation. 
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Further investigation is needed for the following questions: 

• Is A-stability necessary to rule out parametric resonance? 

• What about integration of nonlinear ODEs? 

• How can a "safe" step size oscillation amplitude be estimated for ODE 
integration methods that have numerical damping? 
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